In [1]:
library(Seurat)
library(ggpubr)
library(ggsignif)
library(ggplot2)
obj.integrated=readRDS(file = "../05.2-subtype/scRNA_annodata.KC.rds")
#obj.integrated=subset(obj.integrated,idents = c('SC','BC','SC_BC'))
Idents(obj.integrated) <- "celltype2" #######也就是合并的celltype
obj.integrated=subset(obj.integrated,idents = c('GC_SC','SC','SC_BC','BC'))
obj.integrated2=subset(obj.integrated,idents = c('BC'))
obj.integrated$celltype.g <- paste(Idents(obj.integrated), obj.integrated$group, sep = "_")
Idents(obj.integrated) <- "celltype.g"
table(Idents(obj.integrated))

The legacy packages maptools, rgdal, and rgeos, underpinning this package
will retire shortly. Please refer to R-spatial evolution reports on
https://r-spatial.org/r/2023/05/15/evolution4.html for details.
This package is now running under evolution status 0 

Attaching SeuratObject

Loading required package: ggplot2


    SC_young     BC_young  GC_SC_young  SC_BC_young GC_SC_middle    SC_middle 
        3196         2587          386         1436          320         2373 
   BC_middle SC_BC_middle      SC_aged      BC_aged   SC_BC_aged   GC_SC_aged 
        2803         1312         5276         3928         2457          859
In [2]:
gene=readRDS('bm.pathway.rds')
head(gene,2)
names(gene)
$`basement membrane`
  1. 'ACHE'
  2. 'ACAN'
  3. 'ANG'
  4. 'ANXA2'
  5. 'ANXA2P2'
  6. 'APLP1'
  7. 'DST'
  8. 'ENTPD2'
  9. 'CD151'
  10. 'COL2A1'
  11. 'COL4A1'
  12. 'COL4A2'
  13. 'COL4A3'
  14. 'COL4A4'
  15. 'COL4A5'
  16. 'COL4A6'
  17. 'COL5A1'
  18. 'COL7A1'
  19. 'COL8A1'
  20. 'COL8A2'
  21. 'COL9A1'
  22. 'COL9A3'
  23. 'COL15A1'
  24. 'COL17A1'
  25. 'CST3'
  26. 'DAG1'
  27. 'DLG1'
  28. 'EFNA5'
  29. 'MEGF9'
  30. 'FBLN1'
  31. 'FBN1'
  32. 'FGF9'
  33. 'FN1'
  34. 'HSPG2'
  35. 'TNC'
  36. 'ITGA6'
  37. 'LAD1'
  38. 'LAMA2'
  39. 'LAMA3'
  40. 'LAMA4'
  41. 'LAMA5'
  42. 'LAMB1'
  43. 'LAMB2'
  44. 'LAMB3'
  45. 'LAMC1'
  46. 'LAMC2'
  47. 'LOXL1'
  48. 'LOXL2'
  49. 'MATN2'
  50. 'NID1'
  51. 'NTN3'
  52. 'SERPINF1'
  53. 'PTN'
  54. 'SPARC'
  55. 'TGFBI'
  56. 'THBS2'
  57. 'THBS4'
  58. 'TIMP1'
  59. 'TIMP3'
  60. 'USH2A'
  61. 'VTN'
  62. 'COLQ'
  63. 'ATRN'
  64. 'CASK'
  65. 'TMEFF1'
  66. 'NTN1'
  67. 'ADAMTS1'
  68. 'LAMC3'
  69. 'NID2'
  70. 'LAMB4'
  71. 'NTNG1'
  72. 'TMEFF2'
  73. 'EGFL6'
  74. 'ATRNL1'
  75. 'TINAG'
  76. 'EFEMP2'
  77. 'P3H2'
  78. 'ERBIN'
  79. 'NTN4'
  80. 'SMOC1'
  81. 'SMOC2'
  82. 'VWA1'
  83. 'MMRN2'
  84. 'FRAS1'
  85. 'COL18A1'
  86. 'HMCN1'
  87. 'NTNG2'
  88. 'NTN5'
  89. 'OTOL1'
  90. 'EGFLAM'
  91. 'CCDC80'
  92. 'FREM1'
  93. 'FREM3'
  94. 'NPNT'
  95. 'HMCN2'
  96. 'LAMA1'
  97. 'RELL2'
  98. 'COL28A1'
  99. 'VWA2'
  100. 'FREM2'
  101. 'VWC2'
  102. 'AGRN'
  103. 'AMTN'
$`basement membrane disassembly`
  1. 'CMA1'
  2. 'CTSS'
  3. 'HPN'
  1. 'basement membrane'
  2. 'basement membrane disassembly'
  3. 'basement membrane organization'
  4. 'basement membrane assembly'
  5. 'basement membrane collagen trimer'
  6. 'regulation of basement membrane organization'
In [3]:
nfkb=c('IL1A','TNF','IL1R1','TNFRSF1A','TNFRSF1B','IKBKG','IKBKG','RELA','IKBKG','NFKB1','MYD88','MYD88','MYD88','MYD88','IKBKB',
       'IKBKB','RELA','RELA','TNFAIP3','TNFAIP3','CHUK','IKBKB','MYD88','MAP3K7','IKBKG','TRADD','RIPK1','FADD','MAP3K14','NFKB1',
       'TRAF6','MAP3K1','TAB1','TNFAIP3','NFKBIA','RELA','MAP3K7','MAP3K7','MAP3K7','TRAF6','TAB1')
In [4]:
obj.integrated=AddModuleScore(obj.integrated,features = list(gene[[1]]),name=names(gene)[1])
obj.integrated=AddModuleScore(obj.integrated,features = list(gene[[2]]),name=names(gene)[2])
obj.integrated=AddModuleScore(obj.integrated,features = list(gene[[3]]),name=names(gene)[3])
obj.integrated=AddModuleScore(obj.integrated,features = list(gene[[4]]),name=names(gene)[4])
obj.integrated=AddModuleScore(obj.integrated,features = list(gene[[5]]),name=names(gene)[5])
obj.integrated=AddModuleScore(obj.integrated,features = list(gene[[6]]),name=names(gene)[6])
obj.integrated=AddModuleScore(obj.integrated,features = list(nfkb),name='nfkb')
head(obj.integrated,2)

Warning message:
“The following features are not present in the object: ANXA2P2, NTN3, USH2A, TMEFF1, TINAG, P3H2, ERBIN, NTNG2, OTOL1, FREM3, LAMA1, not searching for symbol synonyms”
Warning message:
“Invalid name supplied, making object name syntactically valid. New object name is basement.membrane1; see ?make.names for more details on syntax validity”
Warning message:
“Invalid name supplied, making object name syntactically valid. New object name is basement.membrane.disassembly1; see ?make.names for more details on syntax validity”
Warning message:
“The following features are not present in the object: NTNG2, not searching for symbol synonyms”
Warning message:
“Invalid name supplied, making object name syntactically valid. New object name is basement.membrane.organization1; see ?make.names for more details on syntax validity”
Warning message:
“The following features are not present in the object: NTNG2, not searching for symbol synonyms”
Warning message:
“Invalid name supplied, making object name syntactically valid. New object name is basement.membrane.assembly1; see ?make.names for more details on syntax validity”
Warning message:
“The following features are not present in the object: OTOL1, not searching for symbol synonyms”
Warning message:
“Invalid name supplied, making object name syntactically valid. New object name is basement.membrane.collagen.trimer1; see ?make.names for more details on syntax validity”
Warning message:
“Invalid name supplied, making object name syntactically valid. New object name is regulation.of.basement.membrane.organization1; see ?make.names for more details on syntax validity”
Warning message:
“The following features are not present in the object: MAP3K14, not searching for symbol synonyms”
A data.frame: 2 × 27
orig.identnCount_RNAnFeature_RNAsamplegroupgroup2percent.mitointegrated_snn_res.1.8seurat_clustersintegrated_snn_res.1.6⋯SC1BC1celltype.gbasement.membrane1basement.membrane.disassembly1basement.membrane.organization1basement.membrane.assembly1basement.membrane.collagen.trimer1regulation.of.basement.membrane.organization1nfkb1
<chr><dbl><int><chr><chr><chr><dbl><chr><fct><chr>⋯<dbl><dbl><chr><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
Skin-Y-18_AAACCTGCACAGTCGC-18Skin-Y-18101283125Skin-Y-18young180.07099131022 ⋯2.192695-0.4098248SC_young-0.03455841-0.028816798-0.1380910-0.03958198-0.01589016 0.05318944-0.13721560
Skin-Y-18_AAACCTGGTCTCAACA-18Skin-Y-18 41451426Skin-Y-18young180.142822689414⋯1.622463 0.6968192SC_young-0.05654325-0.009965443-0.1160495-0.04794820-0.01557398-0.04592599-0.03683826
In [5]:
options(repr.plot.width=12, repr.plot.height=12)
Idents(obj.integrated)=obj.integrated$celltype.g
DefaultAssay(obj.integrated)='RNA'
levels(obj.integrated)=c('GC_SC_young','GC_SC_middle','GC_SC_aged','SC_aged','SC_middle','SC_young','SC_BC_aged',
                         'SC_BC_middle','SC_BC_young','BC_aged','BC_middle','BC_young')
VlnPlot(obj.integrated,features = 'basement.membrane1')
VlnPlot(obj.integrated,features = 'basement.membrane.disassembly1')
VlnPlot(obj.integrated,features = 'basement.membrane.organization1')
VlnPlot(obj.integrated,features = 'basement.membrane.assembly1')
VlnPlot(obj.integrated,features = 'basement.membrane.collagen.trimer1')
VlnPlot(obj.integrated,features = 'regulation.of.basement.membrane.organization1')
VlnPlot(obj.integrated,features = 'nfkb1')
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [6]:
options(repr.plot.width=12, repr.plot.height=12)
Idents(obj.integrated)=obj.integrated$celltype2
DefaultAssay(obj.integrated)='RNA'
levels(obj.integrated)=c('GC_SC','SC','SC_BC','BC')
VlnPlot(obj.integrated,features = 'basement.membrane1')
VlnPlot(obj.integrated,features = 'basement.membrane.disassembly1')
VlnPlot(obj.integrated,features = 'basement.membrane.organization1')
VlnPlot(obj.integrated,features = 'basement.membrane.assembly1')
VlnPlot(obj.integrated,features = 'basement.membrane.collagen.trimer1')
VlnPlot(obj.integrated,features = 'regulation.of.basement.membrane.organization1')
VlnPlot(obj.integrated,features = 'nfkb1')
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

提取重做¶

In [7]:
head(obj.integrated,2)
gene1=c('basement.membrane1','basement.membrane.organization1','basement.membrane.assembly1',
        'regulation.of.basement.membrane.organization1','nfkb1','EXTRACELLULAR_MATRIX_STRUCTURAL_CONSTITUENT_CONFERRING_COMPRESSION_RESISTANCE1',
       'EXTRACELLULAR_MATRIX_CONSTITUENT_CONFERRING_ELASTICITY1','EXTRACELLULAR_MATRIX_STRUCTURAL_CONSTITUENT_CONFERRING_TENSILE_STRENGTH1')
exprs <- data.frame(FetchData(object = obj.integrated, vars = c(gene1,'celltype2','celltype.g')))
A data.frame: 2 × 27
orig.identnCount_RNAnFeature_RNAsamplegroupgroup2percent.mitointegrated_snn_res.1.8seurat_clustersintegrated_snn_res.1.6⋯SC1BC1celltype.gbasement.membrane1basement.membrane.disassembly1basement.membrane.organization1basement.membrane.assembly1basement.membrane.collagen.trimer1regulation.of.basement.membrane.organization1nfkb1
<chr><dbl><int><chr><chr><chr><dbl><chr><fct><chr>⋯<dbl><dbl><chr><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
Skin-Y-18_AAACCTGCACAGTCGC-18Skin-Y-18101283125Skin-Y-18young180.07099131022 ⋯2.192695-0.4098248SC_young-0.03455841-0.028816798-0.1380910-0.03958198-0.01589016 0.05318944-0.13721560
Skin-Y-18_AAACCTGGTCTCAACA-18Skin-Y-18 41451426Skin-Y-18young180.142822689414⋯1.622463 0.6968192SC_young-0.05654325-0.009965443-0.1160495-0.04794820-0.01557398-0.04592599-0.03683826
 
Warning message in FetchData.Seurat(object = obj.integrated, vars = c(gene1, "celltype2", :
“The following requested variables were not found: EXTRACELLULAR_MATRIX_STRUCTURAL_CONSTITUENT_CONFERRING_COMPRESSION_RESISTANCE1, EXTRACELLULAR_MATRIX_CONSTITUENT_CONFERRING_ELASTICITY1, EXTRACELLULAR_MATRIX_STRUCTURAL_CONSTITUENT_CONFERRING_TENSILE_STRENGTH1”
In [18]:
saveRDS(obj.integrated,'scRNA_annodata.KC.BMscore.rds')
In [17]:
noise <- rnorm(n = length(x = exprs[,gene1])) / 100000
exprs[,gene1] <- exprs[, gene1] + noise
#exprs$celltype.g=factor(exprs$celltype.g,levels = c('SC_aged','SC_middle','SC_young','SC_BC_aged','SC_BC_middle','SC_BC_young','BC_aged','BC_middle','BC_young'))
exprs$celltype.g=factor(exprs$celltype.g,levels = c('GC_SC_young','GC_SC_middle','GC_SC_aged','SC_young','SC_middle','SC_aged',
                                                    'SC_BC_young','SC_BC_middle','SC_BC_aged',
                                                  'BC_young','BC_middle','BC_aged'))
#gene1=c('basement.membrane1','basement.membrane.disassembly1','basement.membrane.organization1','basement.membrane.assembly1',
#       'basement.membrane.collagen.trimer1','regulation.of.basement.membrane.organization1')
if(!dir.exists('aging')){
    dir.create('aging')
}
for(i in gene1){
    #compar=list(c('BC_aged','BC_middle'),c('BC_middle','BC_young'),c('BC_aged','BC_young'))
    compar=list(c('BC_aged','BC_middle'),c('BC_aged','BC_young'),c('BC_middle','BC_young'))
    tmp=exprs[,c('celltype.g',i)] ########取出这个通路score，防止ggplot里y输入是字符导致无法计算
    colnames(tmp)=c('celltype.g','pathway')
    plot=ggplot(data = tmp,mapping = aes(x = celltype.g,y = pathway)) +
    geom_violin(scale = "width",adjust =1,trim = TRUE,mapping = aes(fill = celltype.g)) +
    geom_boxplot(width=.2,col="black",fill="white")+
    theme(panel.background = element_blank(),axis.title.x = element_blank(),axis.title.y = element_blank(),
       axis.text.x = element_text(size=10),axis.text.y = element_text(size=10),axis.line = element_line(colour = "black"))+ggtitle(i)+
    stat_compare_means(comparisons = compar,label = "p.signif",method = 't.test',size=10,p.adjust.methods='bonferroni') #,method = 't.test'
    ggsave(paste0('aging/',i,'.pdf'),plot,width = 10,height = 10)
    ggsave(paste0('aging/',i,'.png'),plot,width = 10,height = 10)
}
In [13]:
exprs$celltype.g=factor(exprs$celltype.g,levels = c('GC_SC_young','GC_SC_middle','GC_SC_aged','SC_young','SC_middle','SC_aged',
                                                    'SC_BC_young','SC_BC_middle','SC_BC_aged',
                                                  'BC_young','BC_middle','BC_aged'))
compar=list(c('BC_aged','BC_middle'),c('BC_aged','BC_young'),c('BC_middle','BC_young'))
    tmp=exprs[,c('celltype.g','basement.membrane1')] ########取出这个通路score，防止ggplot里y输入是字符导致无法计算
    colnames(tmp)=c('celltype.g','pathway')
    plot=ggplot(data = tmp,mapping = aes(x = celltype.g,y = pathway)) +
    geom_violin(scale = "width",adjust =1,trim = TRUE,mapping = aes(fill = celltype.g)) +
    geom_boxplot(width=.2,col="black",fill="white")+
    theme(panel.background = element_blank(),axis.title.x = element_blank(),axis.title.y = element_blank(),
       axis.text.x = element_text(size=10),axis.text.y = element_text(size=10),axis.line = element_line(colour = "black"))+ggtitle('basement.membrane')+
    stat_compare_means(comparisons = compar,label = "p.signif",method = 't.test',size=10,p.adjust.methods='bonferroni') 
plot
No description has been provided for this image
In [ ]: